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Abstract 

Mass action type deterministic kinetic models of ion channels 
are usually constructed in such a way as to obey the principle of 
detailed balance (or, microscopic reversibility) for two reasons: first, 
the authors aspire to have models harmonizing with thermodynamics, 
second, the conditions to ensure detailed balance reduce the number 
of reaction rate coefficients to be measured. We investigate a series 
of ion channel models which are asserted to obey detailed balance, 
however, these models violate mass conservation and in their case only 
the necessary conditions (the so-called circuit conditions) are taken 
into account. We show that ion channel models have a very specific 
structure which makes the consequences true in spite of the imprecise 
arguments. First, we transform the models into mass conserving 
ones, second, we show that the full set of conditions ensuring detailed 
balance (formulated by Feinberg) leads to the same relations for the 
reaction rate constants in these special cases, both for the original 
models and the transformed ones. 
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1 Introduction 



1.1 Detailed balancing or microscopic reversibility 

At the beginning of the 20th century it was Wegscheider [27] who gave the 

fcl A ^'^ A 

formal kinetic example A ^ B 2A ^ A + B to show that in some cases 

the existence of a positive stationary state alone does not imply the equality 
of all the individual forward and backward reaction rates in equilibrium: a 
relation = 7^) should hold between the reaction rate coefficients to 

fc-2 ' 

ensure this. Equalities of this kind will be called (and later exactly defined) 
as spanning forest conditions below. Let us emphasize that violation of 
this equality does not exclude the existence of a positive stationary state; it 
exists and it is unique for all values of the reaction rate coefficients, see the 
details in subsection 2.4. 




(a) The Wegscheider reaction (b) The triangle reaction 



Figure 1 

A similar statement holds for the reversible triangle reaction in Fig. lb. 
The necessary and sufficient condition for the existence of such a positive 
stationary state for which all the reaction steps have the same rate in the 
forward and backward direction is kik2k^ = A;_ifc_2^-3- Equalities of this 
kind will be called (and later exactly defined) as circuit conditions below. 
Again, violation of this equality does not exclude the existence of a positive 
stationary state; it exists and is unique for all values of the reaction rate 
coefficients, see the details in subsection 2.4. 

A quarter of a century after Wegscheider the authors Fowler and Milne 
[11] formulated in a very vague form a general principle called the principle 
of detailed balance stating that in real thermodynamic equilibrium all 
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the subprocesses (whatever they mean) should be in dynamic equihbrium 
separately in such a way that they do not stop but they proceed with tlie same 
velocity in both directions. Obviously, this also means that time is reversible 
at equilibrium, that is why this property may also be called microscopic 
reversibility. 

A relatively complete summary of the early developments was given by 
Tolman [23]. 

The modern formulation of the principle accepted by lUPAC [12] 
essentially means the same: "The principle of microscopic reversibility at 
equilibrium states that, in a system at equilibrium, any molecular process 
and the reverse of that process occur, on the average, at the same rate." 

Neither the above document nor the present authors assert that the 
principle should hold without any further assumptions; for us it is an 
important hypothesis the fulfilment of which should be checked individually 
in different models. 

It turned out that in the case of chemical reactions this general principle 
can only hold if both the spanning tree conditions and the circuit conditions 
are fulfilled. However, it became a general belief among people dealing with 
reaction kinetics that the circuit conditions alone are not only necessary but 
also sufficient for all kinds of reactions: Wegscheider's example proving the 
contrary was not known well enough. Vlad and Ross [25] draw the conclusions 
from the Wegscheider example in full generality, but it was Feinberg [10] who 
gave the definitive solution of the problem in the area of formal kinetics: 
he clearly formulated, proved and applied the two easy-to-deal-with sets 
of conditions which together make up a necessary and sufficient condition 
of detailed balance (for the case of mass action kinetics). In other words, 
he completed the known necessary condition (the circuit conditions) with 
another condition (the spanning forest conditions) making this sufficient, as 
well. 

The reason why the false belief is widespread is that in case of reactions 
with deficiency zero the circuit conditions alone are also sufficient not only 
necessary, and most textbook examples have deficiency zero. 

1.2 Ion channel models 

Recent papers on formal kinetic models of ion channel gating show that 
people in this field think that the principle of detailed balance or microscopic 
reversibility should hold. (However, some authors do not consider the 
principle of microscopic reversibility indispensable, e. g. Naundorf et 
al. [19, Supplementary Notes 2, Fig. 3SI(a), page 4] provides a channel 
model which is not even reversible, let alone detailed balanced.) This may 
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be supported either by a theoretical argument: they should obey the 
laws of thermodynamics, or by a practical one: if the principle holds 

one should measure fewer reaction rate coefficients because one also has 
the constraints implied by the principle. The second argument seems to 
be the more important one in the papers by Colquhoun et al. [5], [7]. 
However, the principle is applied in an imprecise way: first, only the 
necessary part consisting of the circuit conditions is applied, second, the 
models are formulated in a way that they do not obey the principle of 
mass conservation. In the present paper we transform the models into 
mass conserving ones, and apply the full set of necessary and sufficient 
conditions. Our main result is that in classes of models including all the 
known ion channel examples are compartmental models, therefore they 
have zero deficiency at the beginning, and being transformed into a mass 
conserving model they have no circuits, therefore one has only to test the 
spanning forest conditions. It is not less interesting that the spanning forest 
conditions obtained for the transformed models are hterally the same as the 
circuit conditions for the original models. 

1.3 Stochastic models 

So far we had in mind only deterministic models (surely not speaking of the 
general but vague formulation of Fowler and Milne). Turning to stochastic 
models one possible approach is to check the fulfilment of microscopic 
reversibility in the following way. Let us suppose we have some measurements 
on a process, and present the data with reversed time, finally use a statistical 
test to see if there is any difference. This is an absolutely correct approach 
and has also been used in the field of channel modeling [21]. 

1.4 Outline 

The structure of our paper is as follows. Section 2 gives a short summary 
of the definitions used and presents Feinberg's theorem. In Section 3 some 
usual ion channel models are transformed into realistic models with mass 
conservation and with the help of a lemma it is shown that in these special 
Ccises the circuit conditions for the origial systems and the spanning forest 
conditions for the transformed systems lead to exactly the same requirements. 
The question of the number of free parameters is also discussed here. Finally 
an outlook and discussion follows in Section 4. The formal proof of our main 
result has been relegated to an Appendix. 

Let us also mention that parts of our investigations has been presented 
in a short, nonrigorous form in [17]. 
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2 Tools to be used 



2.1 Ion channels 

There is a difference in electric potential between the interior of cells and 
the interstitial liquid. An essential part of the system controlling the size of 
this potential difference is the system of ion channels: pores made up from 
proteins in the membranes through which different ions may be transported 
via active and passive transport thereby changing the potential difference in 
an appropriate way. The models of these ion channels are usually described 
in terms of formal reaction kinetics, thus we have to present these notions 
first, then we shall be in the position to present a few alternative models of 
ion channels. 

2.2 Basic definitions of formal kinetics 

Let us consider the reversible reaction 

M M 

^a(m,p)X(m)^5^/?(m,p)X(m) (p = 1, 2, . . . , P) (1) 

m=l m=l 

with M e N chemical species: X(2), . . . , X(M); P e N pairs 

of reaction steps, a{m,p), /3{m,p) e No (m = 1,2, ...,M;p = 
1,2, ...,P) stoichiometric coefficients or moleculairities, and suppose 
its deterministic model 

p 

c'mit) = fm{c{t)) := 5^(/3(m,p) - a{m,p)){w+p{c{t)) - w_p{c{t))) (2) 

p=i 

CmW = c„oeM+ (m = 1,2,...,M) (3) 
— describing the time evolution of the concentration vs. time functions 

t^Cm(t) := [X(m)]{t) 
of the species — is based on mass action type kinetics: 

M 

^+,(c) := fc+,c"(-'^):=A;+,n^?'^''^ (4) 

M 

w.pic) := fc_^c'3(-'^'):=A;_,n<^'^''^ (5) 
{p = 1,2,. ..,F). 
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((2) is also called the induced kinetic differential equation of the reaction 
(1).) The number of complexes is the number of different complex vectors 
among a{.,p) and f3{.,p), i.e. it is the cardinality of the set 

{a{.,py,p = 1, 2, . . . , P} U {/?(., p); p = 1, 2, . . . , P} 

and it is denoted by A^. The Feinberg— Horn— Jackson graph (or, FHJ- 
graph, for short) of the reaction is obtained if one writes down all the complex 
vectors (or simply the complexes, the formal linear combinations on both 
sides of (1)) exactly once and connects two complexes with an edge (or 
two different edges pointing into opposite directions) if there is a reaction 
step taking place between them. Let us denote the number of connected 
components of this graph by L. 

The stoichiometric space is the hnear subspace of generated by 
the reaction vectors: span{^(.,p) — a(.,p);p = 1, 2, . . . , P}; its dimension 
is denoted by S. Finally, the nonnegative integer 5 :— N — L — S is the 
deficiency of the reaction (1). 

Examples to show the meaning of the definitions follow. 

Example 1 (Simple bimolecular reaction) In the simple reversible 
bimolecular reaction A + B ^ C we have Af = 3, P = 1; X{1) = A, X(2) = 
B,X{3) = C; and the complexes are A + B and C, thus the corresponding 
complex vectors are (1,1,0) and (0,0,1). As = 2, L = 1,5' = 1; the 
deficiency of the reaction is 0. 

Example 2 (Triangle reaction) In the triangle reaction (Fig. lb) we have 
M = 3,P = 3;X(1) = A,X{2) = B,X{3) = C; and the complexes are 
A, B and C, thus the corresponding complex vectors are (1, 0, 0), (0, 1, 0) and 
(0, 0,1). As N — 3, L — 1, S — 2; the deficiency of the reaction is 0. 

Example 3 (Wegscheider) In the Wegscheider reaction (Fig. la) we have 
M = 2,P = 2;X{1) = A,X{2) = B; and the complexes are A,B,2A and 
A + B, thus the corresponding complex vectors are (1, 0), (0, 1), (2, 0) and 
(1, 1); therefore the reaction vectors are (1, —1) and (—1, 1). As N — A, L — 
2,S — 1; the deficiency of the reaction is 1. 

Let us mention here that it is a boring task with many possibilities of mistake 
to calculate the characteristic quantities of reactions and this is one of the 
reasons why a program package ReactionKinetics .m is being developed in 
Mathematica, see [24]. The second example may be prepared for the present 
purposes as follows. 
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In[l] 
In [2] 
In [3] 



Out [3] 



< < ReactionKinetics ' 
triangle = {A^B^Cf^A}; 
Column [ReactionsData [triangle] ] 
species— i'{A,B,C} 

ajt->3 

externalspecies— )■{ } 
(£ ^0 

complexes— )-{A , B , C} 

reactionsteps-)- {A-)-B, B-)-A, B-)-C, C-)-B, C-)-A, A-)-C} 
UK ^6 



variables— 7- { c^ , c^ , c c-} 
1 1 
110 
1 1 
10 10 
/3-^[ 10 10 
1 1^ 
-1 1 o' 1 
7^1 1-1-1 1 
^0 1-1-1 
In [4] := ShowFHJGraph [triangle, {ki, k_i, k2, k_2j ks, k_3 }, 
VertexLabeling — > True, DirectedEdges — > True]] 
Other uses of the package are described in the work mentioned above. 



2.3 Models of ion channels 

In the models of ion channels the relevant species are receptors and molecules 
modifying the operation of receptors so as to change the sizes of the pores, 
thereby decreasing or increasing the quantity of ions flowing through the 
channels. Altogether there are several hundreds of different types of ion 
channels in living cells. 

One possible model, see Fig. 2, contains receptors, transmitters, and 
receptor transmitter complexes each with a different conformation having 
different ion- conductance, and these conformations correspond to states in 
which the channels are between the open and closed states [9]. Another 



R+4T T,R' t,r^ z^^^i^ t,r' 

Figure 2: The Erdi-Ropolyi model with four transmitters and three different 
states of the transmitter-receptor complex 
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approach, see Fig. 3, might involve multiple types of modifying molecules 
and complexes, again representing different states of the channels [8] . These 


















■^110 




■^OJO 




"^101 




■^001 




■^111 




'^011 









Figure 3: A model by De Young and Keizer 



are the models we are especially interested in. 

Ion channel models are usually required to fulfil the principle of 
detailed balance both from theoretical and practical points of view. First, 
thermodynamics is said to require the principle to hold, second, if this 
principle holds then the number of reaction rate constants to be measured 
are reduced. Let us turn to the formal definition of detailed balance in the 
framework given in subsection 2.2. 



2.4 Detailed balance: definition and the naive approach 

Within the model exactly defined above we can formulate the property of 
being detailed balanced [14]. Consider the reaction (1) endowed with mass 
action kinetics. 

Definition 1 If G (]R+)*^ is such that 

V*"^"^^ = ^-pc/(-'^^ (p=l,2,...,P), (6) 

then reaction (1) is said to be detailed balanced at the stationary point 

c*. If the reaction is detailed balanced at all its positive stationary points, 
then it is detailed balanced. 

We are especially interested in reactions which are detailed balanced for some 
choices of the reaction rate constants, and also in the restrictions upon the 
rate constants which ensure detailed balancing. 



10 



Example 4 (Simple bimoleculcir reaction) The deterministic model of 

ki 

the reaction A + B C according to subsection 2.2 can be seen to be (in 

k-i 

accord with the usual formulation) 

a' = —kiab + k-ic b' = —kiob + k-ic d = kiab — k-ic 
a(0) = ao 6(0) = bo c(0) = cq 

which simplifies to 

a'{t) — —kia{t){a{t) — ao + bo) + k-i{—a{t) + ao + Co) 

= —kia{t)^ + {kiao - kibo - /c_i)a(t) + k_i{ao + Cq) 

= -k-i{Ka{t)^ -{K{ao-bo)-l)a{t)-ao-CQ) . (7) 

If the reaction starts from nonnegative initial concentrations ao,bo,Co for 
which ao + Co > 0, the unique positive (relatively asymptotically stable) 
equilibrium concentration 

a* = ^{-l + K{ao-bo)+r) 

1 + K{ao + bo + 2co) - r 
* " K{-l + K{ao-bo)+r) 

= -^(l + i^(ao + 6o + 2co) -r) 
ZK 

where K := r := ^/l + 2K{ao + 6o + 2co) + K\ao - bof 

k-i 

will be attained. The reaction is detailed balanced at this vector of stationary 
concentrations for all values of the reaction rate coefficients, i. e. kia^b^ — 
A;_ic* always holds. 

Example 5 (Triangle reaction) The induced kinetic differential equation 
of the reversible triangle reaction being 

a' — —kia + k^ib — k^a + k^c 
b' — kia — k^ib — k2b + /c_2C 
c' — k2b — k-2C + k-^a — k^c 

together with the mass conservation relation 

a{t) + b{t) + c{t) ^ ao + bo + bo ^: m 
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imply that the unique, relatively asymptotically stable vector of positive 
stationary concentrations — if at least one of the initial concentrations 
o-o^bo, Co is positive — are as follows. 

171 

a* = {k_2k-i + ik_i + k2)h)— (8) 

Tn 

K = {k_3k_2 + {k-2 + k3)ki)— (9) 

Tn 

= {k-3k-i + {k-3 + ki)k2)— (10) 

with d := k-2{k-i + ki) + kik2 + k^3{k^2 + k^i + k2) + A;-iA;3 + kik^ + fc2/i;3- 

The reaction is detailed balanced at this vector of stationary concentrations — 
i. e. 

/ciO* — k—ib^ A;2^* — k—2C* k^c^, = k—^a^, 

— if and only if 

kik2k3 — k-ik-2k-3 (H) 

holds. 

Example 6 (Wegscheider) The induced kinetic differential equation of 
the Wegscheider reaction being 

a' — —kia + k-ib — k2a^ + A;_20& 
b' — kia — k-ib + k2a^ — k^2Ci'b 

— which simplifies to 

a' — —kia + k^i{aQ + bQ — a) — k2a^ + A;_2a.(ao + bQ ~ a) 

= -{k2 + k_2)a^ - {ki + k^i - k^2{ao + bo))a + A;_i(ao + bo). 

— together with the mass conservation relation 

a{t) + b{t) = ao + bo =: m 

imply that — unless all the initial concentrations are zero — the unique 
positive (relatively asymptotically stable) stationary concentration vector is 
as follows. 

k_i + ki- k_2m - r ^^^^ 



-2{k_2 + k2) 
+ fci + A;_2?Ti + 2k2m — r 



(13) 



2{k_2 + k2) 

with r := \/{k_i + ki - /c_2m)2 + Ak_im{k_2 + ^2)- (14) 
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The reaction is detailed balanced at this vector of stationary concentrations- 
i. e. 

— if and only if 



holds. 



2.5 The necessary and sufficient condition of detailed 
balancing 

The necessary and sufficient conditions are formulated in the following way 
in [10]. Consider the reaction (1) endowed with mass action kinetics. 

First suppose that we have chosen an arbitrary spanning forest for the 
FHJ-graph of the network. It is possible to find a set oi P — {N — L) 
independent circuits induced by the choice of the spanning forest. For each 
of these circuits we write an equation which asserts that the product of the 
rate constants in the clockwise direction and the counterclockwise direction 
is equal. Thus we have P — [N — L) equations: the circuit conditions. 

Next, these equations are supplemented with the 5 spanning forest 
conditions as follows. Suppose that the edges of the spanning forest has 
been given an orientation. Then there are 5 independent nontrivial solutions 
to the vector equation '^ij^ij = where the sum is taken for all reaction 
steps in the oriented spanning forest and Vjj is the corresponding reaction 
step vector. With these aij coefficients the spanning forest conditions are 

U'^r = U^r, (16) 

where kij are the corresponding rate coefficients. 

With all these the widely- accepted necessciry conditions (the circuit 
conditions) are complemented with the spanning forest conditions to form 
a set of necessary and sufRcient conditions for detailed balancing in mass 
action systems of arbitrary complexity. 

Theorem 1 (Feinberg) The reaction (1) is detailed balanced for all those 
choices of the reaction rate constants which satisfy the P — {N — L) circuit 
conditions and the 6 spanning forest conditions. 

Remark 1 The circuit conditions are called spanning tree method in [7]. 
Remark 2 There are three interesting special cases. 
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1. For a reversible mass action system which has a deficiency of zero, the 
circuit conditions alone become necessary and sufficient for detailed 
balancing. The reason why the circuit conditions were generally 
accepted as sufficient as well, is that a large majority of models are 
of zero deficiency. This case is exemplified by the triangle reaction. 

2. For networks with no nontrivial circuits, that is, in which there are just 
N — L reaction pairs and so P — {N — L) =0, the circuit conditions are 
vacuous. Therefore, the spanning forest conditions alone are necessary 
and sufficient for detailed balancing. The example by Wegscheider 
belongs to this category. 

3. Finally, if a reversible network is circuitless and has a deficiency of 
zero, both the circuit conditions and the spanning forest conditions 
are vacuous. The system is detailed balanced (or fulfils the principle of 
microscopic reversibility), regardless of the values of the rate constants. 
Such is a compartmental system with no circles in the FHJ-graph, the 
simple bimolecular reaction or the Erdi-Ropolyi model. 

3 The main result 
3.1 Our strategy 

Let us denote by M, P, 6, N, L, S, K and M', P', 5', N', L', S', K' the number 
of species, the number of (half) reaction steps, the deficiency, the number of 
complexes, the number of linkage classes, the dimension of the stoichiometric 
space (i.e., the number of independent reaction steps) and the number of 
independent cycles respectively in the original and in the transformed system. 

All the investigated original (not mass-conserving) ion channel models are 
formally compartmental systems which means that each complex consists of 
a single species and all species are different. Therefore all these models are 
of deficiency zero. Thus, in order to check detailed balancing it is enough to 
test the circuit conditions, and this is what the authors in [5, 7] do. 

What we propose is to transform these models into a mass- conserving 
model in such a way as to reflect the same physical reality. The transformed 
models have the following properties. 

1. There is no cycle in the transformed system. 

2. S^S' 

3. N' - L' - S' ^S' ^ K 
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4. The circuit conditions in the original system are equivalent to the 
spanning forest conditions in the transformed system. 

This transformation is constructed in the Appendix for a large class of 
systems — those with rectangular grids as FHJ-graphs — containing all the 
special cases we have met up to now. 

3.2 Lemma 

Consider a directed graph whose edges and vertices are the edges and vertices 
of a planar rectangular grid. Suppose that the graph has n vertices and that 
to each vertex j we assign a vector in R""*"^ such that these vertex vectors 
are linearly independent. Let Ci and C2 be vectors in R""*"^ such that they are 
linearly independent of each other and of each y^. Let us denote by e^j the 
directed edge of the graph from vertex i to vertex j and to each Cjj edge let us 
assign the Vj^ = y^ — y^ vector. Let us define the vectors in the following 
way. If eij is directed in the positive or negative direction in relation to the 
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V32+C1 
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Figure 4 

X axis then Ujj := Vj^ — Ci or Ujj := Vj^ + Ci, respectively. Similarly, if 
is directed in the positive or negative direction in relation to the y axis then 
Ujj := Vjj — C2 or VLij :— Vij + C2, respectively. Let us denote by span{vy} 
the subspace generated by the Vij vectors. 

Lemma 1 Under these conditions the following statements hold. 

1. Along each directed circle in the graph, ^ dij^ij = (^ij'^ij = where 
ttij := 1 if the edges of the graph and the circle are directed in the same 
way and a^j := — 1 otherwise. 
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2. The dimension of span{vij} and span{uij} is n — 1. 

Proof 1. Since the Vjj vectors are the differences of the corresponding 
vertex vectors, it is obvious that along a directed circle, the sum of the v^j 
vectors is 0. It is enough to show that the Ci and C2 vectors disappear in 
the sum of the Uij vectors. In order to see this, first assume that along a 
directed circle we change the direction of the Cij edges so that each is directed 
clockwise. In this case it is obvious that the sum of the Ci and C2 vectors is 
zero since the number of the "+Ci" and "+C2" vectors is equal to the number 
of the "— Ci" and "— C2" vectors, respectively. Then, changing the original 
directions back, the sign of the Ci and C2 vectors changes twice and thus they 
will not appear in the sum. 

2. Let us choose a spanning tree in the graph consisting of n — 1 of the Cjj 
edges. Then the corresponding Vjj vectors are linearly independent and since 
the Ci and C2 vectors are independent of them, the corresponding Uij vectors 
are also linearly independent. 

Remark 3 It is trivial that the statements of the lemma remain true if either 
Ci or C2 is the zero vector, or, if the graph contains edges that are not part 
of a circle. 

Remark 4 The statements of the lemma are also true for graphs consisting 
of A;-dimensional grids (A; > 3), see Fig. 7 and Fig. 8a, 8b as an illustration 
for the three-dimensional case. 

3.3 Examples 

In the next three examples, the left side of the figure shows the original 
system and the right side of the figure shows the transformed system with 
an oriented spanning forest. Both systems are reversible, the arrows show a 
direction needed to write down the spanning forest conditions. The choice of 
the numbering of the species as well as the direction of the reaction vectors 
is arbitrary but in both systems they are chosen correspondingly. 

Example 7 The system in Fig. 5 can be found in [5]. The meaning of 
the species is as follows: The core of the system is obviously a rectangle, 
the additional parts do not mean an extra problem as the reader can easily 
verify it. The original system consists of M = 10 species, N — 10 complexes, 
L — 1 linkage class and it contains two circles while the transformed system 
contains one more species. A, there are A^' = 14 complexes, L' = 3 linkage 
classes and it is circuitless. In order to compare these systems easily, in both 
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cases let us number the species in the same way and let A be the last one, 
that is, 

X{1) := R, X{2) := AR, X(3) := A2R, X(10) := A^F*, X{11) := A. 

Let us assign a vector y-j e M}^ to the ith species so that yij = 1 ii i — j and 
yij — ii i j where i,j — 1, ... ,11 and let a := yn. 

The complex vectors in Fig. 5a are yi, y2, . . . , yio and the corresponding 
reaction vectors are V21 = yi - y2, V23 = ys - y2, • • • , V7,io = yio - Jr- The 
dimension of span{v2i, V23, . . . ,^7,10} is = 9. Thus, the deficiency of this 
system isS — N — L — S = 0. It means that the circuit conditions are 
necessary and sufficient for detailed balancing. The circuit conditions along 
circles 2365 and 4367 are 

^23^36^65^52 — ^32 ^25 ^56 
^43^36^67^74 — ^34^47^76^63 

The complexes in Fig. 5b are numbered as 1, 2, 2a, . . . , 10 and the complex 
vectors are y'^ = yi, y^ = ya, y^„ = ya + a, . . . , y'^o = yio- The reaction 
vectors are U21 = V21, U23 = V23 - a, U43 = V43 + a, U52 = V52, U36 = V36, 

U74 = V74, Uqs = V65 + a, U67 = V67 — a, U58 = V58, Ugg = Vgg, U7_io = V7_io. 

The lemma can be applied to this system with Ci = a and C2 = 0. The 
dimension of span{u2i, U23, . . . , Uy^io} is also S' — 9. Thus, the deficiency is 
5' = 14 — 3 — 9 = 2. Since this system is circuitless, there are two equations 
according to the spanning forest conditions that ensure detailed balancing. 
Along the circles '2365' and '3476' in both systems, V23 + V36 + vgs + V52 = 
U23 + U36 + U65 + U52 = and V43 + V36 + ve? + V74 = U43 + U36 + + U74 = 0. 
Since each coefficient of the Uij vectors in the above linear combinations is 1, 

k' h' h' k' — h' k' k' k' 
'^23 '^36 '^65 ^^'52 ~ '^32'^63'*'56'^25 

Jc' h' h' h' — h' h' h' h' 
'^43'^36'*'67'^74 ~ '>'34 ^^63 ^^76 '*'47 
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which are equivalent to the circuit conditions. 

Remark 5 Let us observe that the equivalence of the circuit conditions in 
the original system and the spanning forest conditions in the transformed 
system follows from the first statement of the lemma, that is, along each 
circle the Vj^ vectors and the correspondingly chosen Uij vectors satisfy the 
same hnear equalities. If, say, instead of circle '4367' we choose circle '234765' 

then V23 - V43 - V74 - V67 + Ves + V52 = U23 - U43 - U74 - U67 + U65 + U52 = 0. 

The corresponding circuit condition in Fig. 5a is 

k23k34k4rkrQke5k52 = ^25 ^56 ^67 ^74^43 
and the equivalent equation from the spanning forest condition in Fig. 5b is 

y (y \-^(y \-^(y \-^y y — y (y \-^(y \-^(y \-^y y 

Example 8 The system in Fig. 6 can be found in [18]. Fig. 6a shows 
the original system where there are M = 11 species and Fig. 6b shows the 
transformed system where there are two more species, A and G. Again, let 
us number the species in the same way as in Fig. 5 and let A and G be the 
last two, that is, 

X{1) := R, X{2) ■=RA,..., X{11) := 6*2/?%, X{12) := A, X(13) := G. 







V12 



@ 



RA 



V32 



© 



RA2 



V4I V25 VS3 

® \® 
RG GRA — ^ GRA2 



V47 VgS V69 

rgY^ G2rX^G2RA2-^ 



(10) 



V9,n 



GiR'Ai 

(a) M = 11,N =11,L= 1, 



(la) 
R+A 



® 
R+G 

"41 

(4) (4a) 

RG RG+A- 

RG+G 

@ 

"47 

RG2 RG2+A- 

® ® 



© 



RA+G 

© 

"25 



RA+A « "^^ RA2 



RA2+G 

©■ 



«63 

© 



GRa ' GRA+A—!^ GRA2 



GRA+G 

© 

«85 



"7^ 



GRA2+G 

© 

"69 



G2RA G2RA+A- 



«98 



G2-RA2- 



"9,10 



© 



G2R*A2 

© 



"9,11 



G2RA2 



(b) M' = 13, N' = 22, L' = 8, 
S' = 10, d' = A,K' = 



Figure 6 
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Let us assign a vector e to the ith species so that y^j — I ii i — j 
and — ii i j where — 13. With a := and g := yi3, the 

corresponding reaction vectors of the transformed system are U12 = V12 — a, 

U32 = V32 + a, U41 = V41 + g, U25 = V25 - g, . . . , Ugg = V98 + a, Ug^io = Vg^io, 

U9,ii = V911. Thus, the lemma can be apphed so that ci = a, C2 = g and 
the X and y axes are directed in the '147' and '123' direction, respectively. 

Example 9 The system in Fig. 7 can be found in [6]. Again, let us number 




(a) M = 18, iV = 18, L = 1, 5 = 17, (b) M' = 20, A^' = 36, L' = 6,S' = 17, 

S = 0,K = 13 S' = 13,K' = 



Figure 7 



the species as shown in Fig. 7a, that is 

X{1) := F*, X{2) := AF*, . . . , X{18) := B2R, X(19) := A, X{20) := B. 

and let us assign a vector yj G to the ith species so that y^j = 1 if 
i = j and y^j = if i 7^ j where i,j = 1,...,20. Similarly as in the 
previous two cases, let Vjj and Ujj respectively denote the reaction vectors 
in the original and in the transformed system (these are the differences of 
the corresponding complex vectors) and let a := yig and b := y2o- Then, 
Uij = Vij — a, Uij = Vij — b or Uij = Vij if Uij and Vij correspond to an edge 
in the graph parallel to the '12', '14' or '17' directions, respectively. The 
lemma can be used here with Ci = a, C2 = b and C3 = 0. 
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Example 10 Consider the system in Fig. 8a where there is one receptor 
with three binding sites, and the different states of the sites are denoted 
by Sijk- The next three figures show three possible transformation of this 
system in the foUowing cases. Fig. 8b shows the transformed versions of 
the system in Fig. 8a in the case when there are three different atoms, A, 
B, C, binding to the three sites. Fig. 8c shows the transformed version of 
the system in Fig. 8a in the case when there are two different atoms. A, 
B, binding to the three sites. This is the De Young and Keizer model, and 
again, the transformed system does not contain a circle. In the interesting 
theoretical case when there is only one atom. A, binding to each of the three 
sites, the transformed system contains a circle, this can be seen in Fig. 8c. 
It can be verified easily that the five circuit conditions in the original system 
are equivalent to the five spanning forest conditions in the systems in Fig. 
8b and 8c and are also equivalent to the four spanning forest conditions and 
one circuit condition in the system in Fig. 8d. 




(a) M = 8, N = 8, L = 1, 
S = 7,6 = 0,K = 5 
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(d) M' = 9,N' = 14, L' = 3, 
S' = 7,6' =4,K' = 1 



Figure 8 
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3.4 On the number of free parameters 



We would also like to make some comments on one of the statements in 
Appendix 2 of [7]. According to this, the number of free parameters can be 
determined as follows. Suppose that we have a system with 

• N complexes, 

• R rate coefficients (as parameters) , 

• and C constraints (the sum of the number of the microscopic 
reversibility constraints and the number of arbitrary constraints — 
independent of the microscopic reversibility constraints and of each 
other — to be imposed on some of the rate coefficients). 

The number of free parameters will then be equal to R — g, where g is the 
rank of an C x matrix, A. 

Recall from [10] that a reversible mass action system is detailed balanced 
if and only if the rate constants satisfy the P — {N — L) circuit conditions 
and the 5 spanning forest conditions where the system has P reaction pairs, 
complexes, L linkage classes, and S is the rank of the stoichiometric space, 
and 6 is the deficiency of the network. Using these notations, it can be written 
that the number of unknowns R equals 2P, and the number of independent 
constraints C equals 

Q + {P-{N-L))+6 = Q + P-S, (17) 

where Q denotes the number of (further independent) external constraints to 
be imposed on some of the rate coefficients. In [10], only P — {N — L)) + 5 = 
P — S is considered to be the number of constraints and in [7], the deficiency 
is not taken into account in this sum. Thus, our equation (17) is a common 
generalization of the equations by Feinberg and Colquhoun et al. 

4 Discussion, open problems 

We have provided a method to transform the most common ion channel 
models into a model where mass-conservation is taken into account. Using the 
theorem by Feinberg we have also shown that the heuristic method happens 
to lead to the same results, in spite of the fact that it is based on imprecise 
assumptions. 

All the original models in question have a rectangular grid structure 
with zero deficiency, and all the transformed models have a deficiency equal 
to the number of independent circuits in the original model. To put it 
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another way, the sum of deficiency and the number of independent circuits is 
invariant under our transformation. The natural question arises if the same 

consequences can be drawn with nonzero deficiency (and nonzero number 
of independent circuits, respectively) and what can be said about reactions 
having an FHJ-graph of different structure. 
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6 Appendix: Reactions of rectangular grid 
structure 



Let us consider a special class of reversible compartmental systems with 
species constructed from D G N different atoms, say, G^, G^, . . . , G^, sitting 
on a receptor which will be omitted as it plays no rule in the calculations. Let 
us represent the species Gl^Gl^ . . . G^^ by the vector (xi, X2, . . . , x^)) G N^, 
and suppose (this is the speciality of the system) that we only have the 
following reaction steps in terms of the atomic representation of the 
species: 

(xi, a;2, . . . ,xd) ^ {xi,X2,. . . , xa-u + 1, Xd+i, ...,xd) (18) 
{0<Xd<Pd-hPdeN;d^l,2,...,D). 

This means that the Feinberg-Horn-Jackson graph (FHJ graph) of the 
reaction is a rectangular grid in the first orthant with H^ill'd + 1) vertex. 

Such kind of reactions are often used when modeling ion channels see Fig. 
8 or [8]. 

Realizing that atoms are not conserved in the above reaction, we try to 
improve it by constructing a model without this fault but reflecting the same 
physical reality. In order to do so we have to introduce D new, single-atom 
species, G"^ {d— 1, 2, . . . , D) and the new reaction steps 

Gd + (Xi, 0:2, ... , Xd) ^ {Xu X2, . . . , + 1, . . . , X£i), (19) 

where e^ is the c/th element of the standard base. 

To test if a general reaction is detailed balanced or not one has to 
write down 5 number of circuit conditions and K number of spanning 
forest conditions in terms of the reaction rate constants which form a set 
of necessary and sufficient conditions together. 

If we are interested in detailed balancing of the first reaction (18) 
we should rather transform it to (19) and have only the spanning forest 
conditions. The astonishing fact, however, is that for these special reactions 
not only the number of conditions are the same, but the conditions 
themselves, as well. 

Let us use the following notations: 
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N : the number of complex vectors (the number of vertices) 

P : the number of reaction pairs (the number of edges) 

L : the number of hnkage classes (the number of connected components) 

S : the dimension of the stoichiometric space 

(the number of independent reaction steps) 

5 : the deficiency 

K : the number of independent circuits 

To get some experience with this kind of systems we summarize the 
essential characteristics of these systems in two and three dimensions and 
then formulate and prove the general formula. 

Statement 1 If D = 2 then the formulas for system (18) (left column) and 
system (19) (right column) are the following: 

N = (p+l)(g + l) N' = 2{p + q) + 3pq 

L = 1 L' — p + q + pq 

S ^ N-1 S' ^ S 

S ^ N - L- S = S' = N'-L'-S' = K 

K = P-{N -L)^pq K' = P' - (AT' - L') = 

where P = P' = p(g + 1) + (p + l)q. 
If D = 3 then the formulas are 



N = (p+l)(? + l)(r + l) N' 

L = 1 L' 

S ^ N-l S' 

5=0 5' 

K — pq + pr + qr + 2pqr K' 



2{p + q + r) + 3{pq + pr + qr) + Apqr 

{p + q + r) + {pq + pr + qr) + pqr 

S 

K 





where P = P' = p{q + l){r + 1) + {p + l)q{r + 1) + {p + l){q + l)r, see Fig. 
8a and 8b as an illustration. 

Theorem 2 

1. The essential characteristics of reactions (18) (with its FHJ-graph as a 
rectangular grid) and (19) are as follows. 
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= 
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5' 


= K 
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k=2 


■■Pd, K' 


= 



where each sum is taken with the restrictions 1 < di < d2 < ■ ■ ■ < < 
D. 

2. The circuit conditions for reaction (18) are exactly the same as the 
spanning tree conditions for reaction (19). 

Proof. In both systems the number of edges can be calculated as 

P^P' = P1(P2 + 1)...(PD + 1) + (P1 + 1)P2(P3 + 1)...(PD + !) + ••• + 

+{pi + 1) . . . {pd-1 + 1)pd = 

D 

= XI kPdiPd2 ■■■Pdk {I < di < d2 < ■ ■■ < dk < D) 

k=l 

The number of independent circuits in a graph can be calculated as K — 
P — [N — L). Thus, using that iV = 1 + L', we obtain the formula for K: 

K = P - {N - L) ^ P - N + 1 ^ P - L' 

D D 

= Yl kPdiPd2 ■■■Pdk-Yl P<^^P<^^ ■■■P<^k 

k=l k=l 

The formulas for A^' and L' follow from the following observation: in the 
graph of the transformed system the number of components consisting of 
one edge (and two vertices) is pi + p2 + ■ ■ ■ +Pd', the number of components 
consisting of two edges (and three vertices) is piP2 +P1P3 + ■ ■ -Pd-iPd', etc.; 
the number of components consisting of D edges (and D + 1 vertices) is 
P1P2 ■ ■ -Pd- The equality S = S' and the equivalence of the circuit conditions 
and spanning forest conditions follow from the D dimensional version of the 
lemma. Finally, using that S' = S = N - 1 = L' , 6' = N' - L' - S' and 
K' ^ P' - {N' - L'), we obtain the formulas for S' and K'. 
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